Extreme Mass-Ratio Inspirals in the Effective-One-Body Approach: 
Quasi-Circular, Equatorial Orbits around a Spinning Black Hole 
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We construct effective-one-body waveform models suitable for data analysis with LISA for 
extreme-mass ratio inspirals in quasi-circular, equatorial orbits about a spinning supermassive black 
hole. The accuracy of our model is established through comparisons against frequency-domain, 
Teukolsky-based waveforms in the radiative approximation. The calibration of eight high-order 
post-Newtonian parameters in the energy flux suffices to obtain a phase and fractional amplitude 
agreement of better than 1 radian and 1% respectively over a period between 2 and 6 months de- 
pending on the system considered. This agreement translates into matches higher than 97 % over 
a period between 4 and 9 months, depending on the system. Better agreements can be obtained 
if a larger number of calibration parameters are included. Higher-order mass ratio terms in the 
effective-one-body Hamiltonian and radiation-reaction introduce phase corrections of at most 30 
radians in a one year evolution. These corrections are usually one order of magnitude larger than 
those introduced by the spin of the small object in a one year evolution. These results suggest 
that the effective-one-body approach for extreme mass ratio inspirals is a good compromise between 
accuracy and computational price for LISA data analysis purposes. 
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I. INTRODUCTION 

Extreme mass-ratio inspirals (EMRIs) are one of the 
most promising sources of gravitational waves (GWs) ex- 
pected to be detected with the proposed Laser Interfer- 
ometer Space Antenna (LISA) [l|-|4| . These sources con- 
sist of a small compact object, such as a neutron star or 
stellar-mass black hole (BH), in a close orbit around a 
spinning, supermassive BH |5[. Gravitational radiation 
losses cause the small object to spiral closer to the su- 
permassive BH and eventually merge with it. Hence, the 
GW signal from such events encodes information about 
strong gravity, allowing tests of general relativity and 
of the Kerr metric |7Hl7| , as well as measurements of the 
spins and masses of massive BHs fl8j . 

Unfortunately, EMRIs are very weak sources of GWs 
at their expected distances from us, and thus, they must 
be observed over many cycles to be detectable § . For ex- 
ample, a typical EMRI at a distance of 3 Gpc would pro- 
duce GWs with signal-to-noise ratios (SNRs) on the order 
of 10-200 depending on the observation time. Therefore, 
matched filtering is essential to extract EMRIs from LISA 
noise and the foreground of unresolved GWs from white 
dwarf binaries in our galaxy. 

Matched filtering consists of cross-correlating the data 
stream with a certain noise-weighted waveform tem- 
plate If the latter is similar to a GW event hidden 
in the data, then this cross-correlation filters it out of the 
noise. Of course, for matched filtering to be effective, one 



must construct accurate template filters. Otherwise, real 
events can be missed, or if an event is detected, parameter 
estimation can be strongly biased [20[ ■ The construction 
of accurate EMRI waveforms is extremely difficult due to 
the long duration of the signal and the strong-field nature 
of the orbits. A one- year EMRI signal contains millions 
of radians in phase information. To avoid significant de- 
phasing, its waveform modeling must be accurate to at 
least one part in 10 5 -10 6 [U. 

Such an exquisite accuracy requirement is complicated 
further by the strong- field nature of the orbit. An EMRI 
can reach orbital velocities of two-thirds the speed of 
light and orbital separations as small as a few times 
the mass of the supermassive companion. This auto- 
matically implies that standard, post-Newtonian (PN) 
Taylor-expanded waveforms fail to model such EMRI or- 
bits [111]. PN theory relies on the assumptions that all 
orbital velocities are much smaller than the speed of light 
and that all objects are at separations much larger than 
the total mass of the system (23|. A better approxima- 
tion scheme to model EMRIs is BH perturbation theory, 
where one only assumes that the mass ratio of the sys- 
tem is much less than unity (24[. This is clearly the case 
for EMRIs, as the mass ratio is in the range 10 -4 -10~ 6 . 
Perturbation theory, however, is computationally and an- 
alytically expensive. Only recently have generic orbits 
been computed around a non-spinning BH to linear or- 
der in the mass ratio [25|, (2(| , and it is unlikely that these 
will be directly used for EMRI data analysis [|| . 
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EMRIs involve complicated inspiral analysis, but un- 
like comparable-mass coalescences, the merger and ring- 
down phase can be completely neglected. To see this, 
note that the instantaneous amplitude of the waves from 
a binary scales as /z, where /i = mi m%/M is the reduced 
mass, M is the total mass and mi. 2 are the component 
masses. The inspiral lasts for a time ~ and releases 
an energy flux ~ /x //i ~ /i. In contrast, the merger and 
ringdown last for a time ~ M (independent of /z), and 
thus, release an energy flux ~ /j 2 . For an EMRI, 
and the inspiral clearly dominates the signal. Based on 
this argument, we neglect the merger and ringdown, fo- 
cusing on the inspiral for our analysis. 



A. Summary of Previous Work 

The modeling of EMRIs has been attempted in the 
past with various degrees of success. One approach is 
to compute the self-field of the test particle to under- 
stand how it modifies the orbital trajectory. This task, 
however, is quite involved, both theoretically and com- 
putationally, as the self-field contains a diver gent piece 



that is difficult to regularize (see, e.g. Rcf. [27j for a 
recent review). Recently, a breakthrough was achieved, 
with the full calculation of the self-force for _generic EM- 
RIs about non-spinning supermassive BHs [25|, [26| . Such 
calculations, however, are computationally prohibitive if 
the goal is to populate a waveform template space. 

Another approach is to use more approximate meth- 
ods to model the EMRI trajectories. One such approach 
was developed by Hughes [H, [2t|, following the pio- 
neering work of Poisson [30j . In this radiative- adiabatic 
scheme, the inspiral is treated as a sequence of adiabati- 
cally shrinking geodesies. The degree of shrinkage is de- 
termined by solving the Tcukolsky equation on each indi- 
vidual geodesic. Its solution encodes how the constants of 
the motion (the energy, angular momentum and Carter 
constant) change due to GW emission. By interpolat- 
ing across such sequence of geodesies, one then obtains 
a continuous inspiral and waveform. The calculation of 
a single waveform, however, is rather computationally 
expensive, as it requires the mapping of the entire or- 
bital phase space, which for generic orbits is likely to be 
prohibitive. It is also worth noting that the radiative ap- 
proximation neglects the impact of conservative effects 
which, especially for eccentric orbits, are likely to be im- 
portant [3l| . 

Other, perhaps more rough approximations can also be 
used to model EMRIs. The templates obtained through 
these methods are sometimes called kludge waveforms to 
emphasize their approximate nature. The goal of their 
construction was never to provide sufficiently accurate 
templates for real data analysis. Instead, kludge wave- 
forms were built to carry out descoping or parameter 
estimation studies to determine roughly the accuracy to 
which parameters could be extracted, given an EMRI de- 
tection with LISA. 



The first kludge waveforms were constructed by Barack 
and Cutler [l8[ . These waveforms employ the quadrupole 
formula to build templates as a function of the orbital tra- 
jectories. The latter are simply Keplerian ellipses with 
varying orbital elements. The variation of these is de- 
termined by low-order PN expressions, constructed from 
the GW energy and angular momentum fluxes. An im- 
provement of these fluxes was developed by Gair and 
Glampedakis [32[ , who fitted these low-order PN expres- 
sion to more accurate fluxes constructed from solutions 
to the Teukolsky equation. A further improvement was 
developed by Babak et al. (33|, who modeled the wave- 
forms via a quadrupole-octopole formula and the orbital 
trajectories via solutions to the geodesic equations, aug- 
mented with PN-orbit-averaged evolution equations for 
the orbital elements. 

All of these improvements, however, do not mean that 
kludge waveforms would be effectual or faithful for real- 
istic data analysis with LISA. One cannot exactly quan- 
tify this statement because exact EMRI waveforms are 
not available and will not be in the near future. One 
can nonetheless predict that these approaches will be in- 
sufficient because critical components of the fluxes are 
not being taken into account. For example, GWs do not 
only escape to infinity, but they are also absorbed by 
the supermassive BH, contributing to the overall fluxes 
of energy and angular momentum. This contribution is 
non-negligible if one considers sufficiently long waveforms 
(longer than a few weeks). In fact, as we shall show in 
this paper, even the inclusion of such terms and very high 
order PN expressions in the fluxes is still insufficient for 
accurate waveform models that last more than a couple 
of months. 



B. The Effective-One-Body Approach 

The effective-one-body (EOB) formalism was intro- 
duced in Refs. [34ll35l] to model the inspiral, merger, and 
ringdown of comparable-mass BH binaries. This scheme 
was then extended to higher PN orde rs 13611 . spinning 
BHs [37l - l40j . small mass-ratio mergers 4l]-43[, and im- 
proved by resumming the radiation reaction-force and 
waveforms [HI, |43 - |47| . In the comparable mass case, 
phase and amplitude agreement was achieved between 
EOB and numerical-relativity waveforms, after calibrat- 
ing a few parameters }48| - |5(J |. By calibrating the EOB 
model to the comparable mass case, one can also im- 
prove the agreement of the model with the self- force pre- 
dictions [25|, [5l| . The combination of EOB and BH per- 
turbation theory tools for LISA data-analysis purposes 
was first carried out in Refs. [H2, [53[. In these papers, 
the EOB scheme was found successful for the coherent 
modeling of EMRIs about a non-spinning background for 
a 2 year period. Here we extend these results to non- 
preccssing EMRIs about a spinning background. 

As a first step toward the construction of accu- 
rate EMRI waveforms, we concentrate on quasi-circular, 
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equatorial EMRIs about a spinning, supermassive Kerr 
BHs. The modeling of such EMRIs is simpler than that 
of inclined and eccentric ones, as only a single compo- 
nent of the radiation-reaction force is non- vanishing and 
entirely controlled by the GW energy flux (the Carter 
constant vanishes by symmetry). Moreover, such EMRIs 
are expected in at least one astrophysical scenario [54| . 
In this setup, stellar-mass compact objects are either cre- 
ated in the accretion disk surrounding the supermassive 
BH or are captured by the disk, and hence move with 
the disk. The accretion disk is expected to be in the spin 
equatorial plane within a few hundred gravitational radii 
of the supermassive BH [55| . 

We first calibrate the EOB energy flux to the leading- 
order energy flux computed in BH perturbation theory 
through the solution to the Teukolsky equation. This cal- 
ibration is more complicated than for non-spinning sys- 
tems because it must now be performed globally, i.e., as 
a function of both spin and velocity. This increases the 
computational cost of the calibration and the number of 
calibration parameters, as a bivariate series generically 
contains more terms than a monovariate one. After cal- 
ibrating 8 parameters, we find that the fluxes agree to 
within one part in 10 3 for all spins [a/A/ = (-0.99,0.99)] 
and velocities [v = (0.01 c, Wisco)] considered, where n ISCO 
is the velocity at the innermost circular orbit (ISCO). 

Once the energy flux has been calibrated, we evolve 
the Hamilton equations in the adiabatic approximation 
and compare the amplitude and phase evolution to that 
obtained with an approximate BH perturbation theory, 
numerical result. For the latter, we employ the so-called 
radiative approximation [28l . [23 | , where one models the 
EMRI as an adiabatic sequence of geodesies with vary- 
ing orbital elements, as prescribed by the solution to 
the Teukolsky equation. We find that the EOB and 
Teukolsky-based waveforms agree in phase and relative 
amplitude to better than 1 radian and 1% respectively 
after 2 or 6 months of evolution, depending on the sys- 
tem considered. Better agreements can be obtained if a 
larger number of calibration parameters were included. 

Our EOB waveforms differ from previous kludge mod- 
els on several fronts. First, the radiation-reaction force is 
here computed differently than in the kludge approach. 
In the latter this force is calculated from PN, Taylor- 
expanded fluxes that encode the GW that escape to 
infinity only. These fluxes were then improved by fit- 
ting a very large number of parameters to more accurate 
Teukolsky-fluxes with a log-independent, power-series ex- 
pansion for the fitting functions [HJ. In the EOB ap- 
proach, the radiation-reaction force is computed directly 
from the factorized resummed waveforms p5l |46| . These 
are enhanced through the addition of BH absorption 
terms and then the calibration of eight high PN-order 
parameters to Teukolsky fluxes with a log- dependent, 
power-series expansion for the fitting functions. Second, 
the conservative dynamics are also treated here differ- 
ently than in the kludge approach. In the latter, the 
Hamiltonian is either a two-body, Newtonian one [lH 



or the full test-particle limit one, i.e., Schwarzschild or 
Kerr [33[ . In the EOB approach, the conservative dynam- 
ics not only encodes the exact test-particle limit Hamilto- 
nian, but they also allow for the inclusion of finite mass- 
ratio terms and of the spin of the small body. 



C. Data Analysis Implications 

The waveforms computed here are thus suitable for 
coherent data analysis over periods of several months. 
This can be established by computing the overlap be- 
tween the EOB and Teukolsky-based waveforms, after 
maximizing over extrinsic parameters (an overall phase 
and time shift). We find that, when eight calibration pa- 
rameters are used, the overlap remains higher than 97 % 
over 4 to 9 months of evolution, depending on the sys- 
tem considered. This is to be compared with numerical 
kludge waveforms [33| whose overlap drops to 56 % and 
74 % after 4 and 9 months respectively, even when forty- 
five calibration parameters are used to fit the flux [32J . 
Of course, one could obtain higher overlaps by maximiz- 
ing over intrinsic parameters, such as the chirp mass or 
the spin of the background, but this would naturally bias 
parameter estimation. Also, when integrating over only 
two weeks, the overlap increases, remaining higher than 
0.99999 at 1 Gpc regardless of the model used. 

The benefit of coherently integrating over longer peri- 
ods of time is that the recovered SNR naturally increases, 
thus allowing us to detect signals farther out and improv- 
ing parameter estimation. One can see this by simply 
noting that the SNR scales with the square root of the 
time of observation. For example, coherent integration 
over 4 or 9 months instead of two weeks increases the 
SNR at 1 Gpc from 13 to 64 and from 27 to 85 for two 
prototypical EMRIs. Such a large increase in SNR by 
coherently integrating over long observation times brings 
EMRIs not only to a confidently detectable range, but 
would also allow interesting tests of GR. 

We conclude the paper by studying the error intro- 
duced in these waveforms due to neglecting second-order 
mass-ratio terms in the radiation-reaction force (dissipa- 
tive PN self-force) and first-order in the mass-ratio terms 
in the Hamiltonian (conservative PN self- force). Such 
mass-ratio dependent effects can easily be included in 
the EOB prescription, as they are known in the PN/EOB 
framework. Of course, since these are known to finite PN 
order, we cannot include full second-order effects. These 
effects should be considered estimates, since the complete 
result may differ from the PN prediction. We find that 
such PN radiation-reaction effects modify the phase of 
the waveform by O(10) radians in a one year evolution, 
provided the EMRI samples the strong-gravity regime 
close to the ISCO. In a two-month period, however, the 
inclusion of finite mass-ratio effects increases the mis- 
match from 2.9 x 10~ 5 to 3.6 x 10 -5 at 1 Gpc. This im- 
plies that such effects will only be seen if one coherently 
integrates over a sufficiently long time of observation. We 
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find a somewhat smaller final dephasing when we allow 
the second body to be spinning and neglect any self-force 
corrections. The relative importance of the conservative 
or dissipative PN self-force terms and that of the spin 
of the second body depends on somewhat on the EMRI 
considered. Generically, we find all such corrections to 
be larger than one radian after a full year of integration, 
while they are negligible over a two month period. 

This paper is organized as follows. Section [TTI describes 
how we model EMRIs analytically and numerically. Sec- 
tion IIHI discusses how the analytical EOB model is cal- 
ibrated to the Teukolsky energy flux. Section IIVI com- 
pares EOB evolutions to Teukolsky ones, while Sec. |V] 
discusses the data analysis implications of such a compar- 
ison. Scction lvTl estimates the effect of mass-ratio depen- 
dent effects and Sec. IVIll concludes and points to future 
research. Appendix |A~1 presents details on the transforma- 
tion between spheroidal and spherical tensor harmonics. 
Appendix |B] contains expressions for the GW energy flux 
absorbed by BHs. Finally, in Appendix [C] we write the 
EOB Hamiltonian derived in Ref. [4(| when BHs carry 
spins aligned or anti-aligned with the orbital angular mo- 
mentum. We use geometric units, with G = c = 1, unless 
otherwise noted. 



II. EMRI MODELING 

A. Analytical modeling: EOB-based waveforms 

Consider a BH binary system with masses mi and 7712, 
total mass M = mi + m 2 , reduced mass /1 = mim 2 /M 
and symmetric mass ratio v = fi/M. We assume that 
the orbital angular momentum is co-aligned or counter- 
aligned with the individual BH spins Sa = aAin-A = 
qAm 2 A , where a a = SaI^a denotes the Ath BH's spin 
parameter and qa = a a/tu a denotes the dimensionlcss 
spin parameter. 

We first discuss the case of a non-spinning BH (92 = 0) 
with mass 7772 orbiting a spinning BH with spin param- 
eter qi and mass mi ^> 777,2, to leading order in the 
mass ratio m 2 /mi. Sublcading terms in the mass ra- 
tio and terms proportional to q 2 introduce conservative 
corrections that are not included in the Teukolsky wave- 
forms, which we shall use to calibrate our model, and we 
therefore neglect them during the calibration. Eventu- 
ally, however, we shall turn these conservative terms on 
and estimate their effect using the spin EOB Hamiltonian 
of Ref. [H (see Appendix [C]). 

In the EOB framework, the orbital trajectories are ob- 
tained by solving Hamilton's equations, supplemented by 
a radiation-reaction force describing the backreaction of 
GW emission on the orbital dynamics. Neglecting con- 
servative corrections of order 0(777,2/777,1) and 0(q 2 ), the 
spin EOB Hamiltonian reduces to the Hamiltonian of a 



non-spinning test-particle in Kerr, i?NS : 

HEOB = H m [1 + 0(7772 / 777!) + 0(q 2 )] , (1) 

H NS = P i Pi + a^jml + -yij Pl Pj , (2) 



where 
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g^u being the Kerr metric. In Boyer-Lindquist coordi- 
nates (i, r, 9, <p) and restricting ourselves to the equatorial 
plane = 0, the relevant metric components read 
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where Wfy = 2qi m\ r, and the metric potentials arc 

A = r 2 + q\ m\ — 2mi r , 

A = (r 2 + q\ m\f - qfmfA. 



(7) 
(8) 



Although the EOB formalism includes possible non- 
adiabaticities in the last stages of the inspiral and plunge, 
it is not necessary to include non-adiabatic effects here. 
Generically, for the systems that we consider, we find that 
the inclusion of non-adiabatic corrections leads to small 
phase corrections (of 0(1 rad) after one year of evolu- 
tion) [52|, [53| . The assumption of adiabaticity allows us 
to simplify the evolution equations that are solved nu- 
merically. This in turn reduces the computational cost of 
producing EOB waveforms: an adiabatic EOB evolution 
requires a few CPU seconds, while a non-adiabatic one 
would require CPU days or weeks. The non-adiabatic 
model is computationally more expensive because one 
needs to solve all of Hamilton's equations, with radiation 
reaction source terms that are expensive to evaluate. 

The Hamiltonian of Eq. simplifies drastically when 
we consider circular, equatorial orbits (8 = tt/2) with 
Si co- aligned or counter-aligned with the orbital angular 
momentum (see eg. [56|). Imposing p r = 0, which is 
a necessary condition for circular orbits, and inserting 
Eqs. (|B"a|) - (|6dp in Eq. @, a straightforward calculation 
returns 
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where 



Q = 1 



2 2 

m\ A 



(10) 



and = L is the conjugate momentum to the <j> Boyer- 
Lindquist coordinate or simply the orbital angular mo- 
mentum. Imposing the condition p r = (dH^§/dr) Pr= Q = 
0, which is also satisfied by circular orbits, we can solve 
for I as a function of r and qi [56j 



L = ±1712 TTli 



1/2 



r 2 =F 2gi m\ /2 r 1 / 2 + q\ m\ 
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,3/4 ^ r 3/2 _ 3 mir l/2 ± 2?1 m J 



1/2 



(11) 

where ± corresponds to prograde or retrograde orbits, re- 
spectively. Inserting the above equation in Eq. © yields 
an expression for the energy E = H^g c (r) of circular or- 
bits in Kerr [Ml 



E = mi 



1 



mi 



1 — Zmi/r ± 2gi m^ 2 /r 3 / 2 



(12) 



The above quantities can also be expressed in terms of 
the orbital velocity u>, once r(ui) is derived for circular 
orbits. Computing uj = (dH^/dp^) Pr= Q, using Eq. (jTTJ) , 
we obtain 



[1 - qi (miw)] 2/3 



(miw 



,2/3 



(13) 



We also define the parameter v = (miw) 1 ' 3 . 

In the adiabatic approximation, the orbital evolution 
is fully determined by the frequency evolution through 
Eq. (|T3|) . Assuming the motion follows an adiabatic se- 
quence of quasi-circular orbits, we can use the balance 
equation L = E /uj = —J- /uj to derive 



(14) 



where T is the GW energy flux (see e.g. Ref. [35|). The 
multipolar factorized form of this flux, proposed in the 
non-spinning case in Refs. [42l . |45| and extended to the 
spin case in Ref. [46|, is given by 



^ 8 l 
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which under the assumption of adiabaticity reduces to 
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with 
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where e p denotes the parity of the multipolar waveform 
(i.e., e p = if t + m is even, e p = 1 if £ + m is odd), and 
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When spin effects are present, the expressions for all the 
terms in Eq. ([17]), namely S^jv), Tg m (v), 8t m (v) and 
pt m (v) can be read in Ref. [46| [see Eqs. (24), (25), (26) 
and (29) therein]. The functions Yi tm (8,<p) are the stan- 
dard spherical harmonics, while nf^ and ce+ ep are nu- 
merical coefficients that depend on the mass ratio (see 
Eqs. (5)- (7) in Ref. [45|). As before, we work to lead- 
ing order in v initially, and later study how the terms of 
higher-order in v affect the GW phase evolution. 

The solution to Eq. (|14[) requires that we prescribe ini- 
tial data. We here choose post-circular initial conditions, 
as described in Ref. , to set-up a mock evolution that 
starts at a separation of lOOmi and ends at either the 
ISCO or whenever the GW frequency reaches 0.01 Hz. 
This mock evolution is then used to read initial data one- 
year before the end of the mock evolution. This approach 
leads to an accurate initial data prescription, without any 
eccentricity contamination. For example, the error in the 
initial frequency induced by starting the mock evolution 
at 100 mi, instead of 200 mi, is on the order of 10 -9 Hz, 
which leads to a difference in accumulated GW cycles of 
0.03 rads after a one year evolution. 

Finite mass-ratio corrections can be incorporated 
into the EOB model by including subleading terms of 
0(1712/ mi) and 0(52) in the Hamiltonian, angular mo- 
mentum and r(uj) relation. We shall first ignore such 
terms to compare against Teukolsky-based waveforms. In 
Sec. IVI1 we shall study how our results change when we 
include such terms. To do so, we shall still assume cir- 
cular, equatorial orbits and an adiabatic evolution, but 
employ the spin EOB Hamiltonian of Refs. [4(], [HtJ (re- 
viewed in Appendix [C]), instead of the Kerr Hamiltonian 
of Eq. ®. 

Except for this change, the EOB waveform model- 
ing with finite mass-ratio corrections follows closely the 
derivation presented above. First, we compute the an- 
gular momentum associated with the Hamiltonian of 
Eq. (|C22[) for circular, equatorial orbits, imposing p r — 
(<9i?EOB/<9r) Pi . = o = and solving for L = p^. Then, 
we derive the orbital frequency ui = (dHEOB/dp<f,) Pl , = o 
to relate r tow, and to express L in terms of uj. When 
mass-ratio corrections are present, however, the Hamil- 
tonian becomes much more involved, so solutions for L 
as a function of r must be searched numerically. We have 
checked that the discretization and interpolation used to 
solve these equations numerically do not introduce an 
error larger than 10~ 10 in the Hamiltonian (|C22|) . 
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B. Numerical modeling: Teukolsky-based 
waveforms 

Teukolsky-based waveforms is the name we give to ra- 
diative models that use the Teukolsky equation to pre- 
scribe the radiation-reaction force for an inspiral. We use 
BH perturbation theory, considering a background space- 
time with mass mi and spin |Si| = m\q (recall that the 
spin parameter q\ — ai/m{). The inspiraling object is a 
test-body with mass mi -C rti\ and no spin (52 = 0). In 
principle, the masses and spins used here should be the 
same as those introduced in the EOB model. 

The radiative approximation assumes that EMRIs can 
be modeled as an adiabatic sequence of geodesies with 
slowly-varying constants of the motion. Consider the 
discretization of the orbital phase space, each point of 
which represents a certain geodesic with a given set of 
constants of motion (energy E, angular momentum L 
and Carter constant Q). For quasi-circular, equatorial 
EMRIs, the Carter constant vanishes, while the varia- 
tion of the energy and angular momentum are related 
via SL = (5E)/oj, us being the orbital frequency. At each 
point in the orbital phase space, the geodesic equations 
can be solved to obtain the orbital trajectory of the small 
compact object, given any spin of the background. 

Once the geodesic trajectories are known, one can use 
these to solve the linearized Einstein equations and ob- 
tain the gravitational metric perturbation. This is best 
accomplished by rewriting the linearized Einstein equa- 
tions in terms of the Newman-Penrose curvature scalar 
■04 to yield the Teukolsky equation [58]. One can de- 
compose -04 into spin- weight —2 spheroidal harmonics 
-2<5fm(0)) using Boyer-Lindquist coordinates (t,r,6,<f)), 
in the Fourier domain: 
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The radial functions Ri muj (r) satisfy the radial Teukolsky 
equation 



A -r- 
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where A is given in Eq. (J7)) and the radial potential is 



V(r) ee - 
with K ee (r 2 



K 2 + 4i(r - mi)K 



A 



+ 8iujr + X , 



(21) 



m 2 q 2 )u> — mniiqi, A = Ei m — 2m\q\muj + 
mfqfuj, and £g m the spheroidal harmonic eigenvalue. The 
source function Timu is given explicitly in Eq. (4.26) of 
Ref. [HI and it depends on the stress-energy tensor for a 
test-particle in a geodesic trajectory. 

The Teukolsky equation admits two asymptotic solu- 
tions: one outgoing as r — ¥ 00 and one ingoing as one 
approaches the background's event horizon. These two 



solutions represent outgoing radiation at future null infin- 
ity and ingoing radiation that falls into the BH through 
the event horizon. Both types of radiation are critical 
in the modeling of EMRIs; not including BH absorption 
can lead to errors in the waveform of order 10 4 radi- 
ans 0, [59| . These solutions can then be used to recon- 
struct both the GW radiated out to infinity, as well as the 
total energy flux lost in GWs. The energy flux can then 
be related to the temporal rate of change of the orbital 
elements, such as the orbital radius. 

Solving the Teukolsky equation for a geodesic orbit 
tells us how that orbit tends to evolve due to the dissipa- 
tive action of GW emission. By doing so for each point 
in orbital phase space, we endow this space with a set 
of vectors that indicate how the binary flows from one 
orbit to another. We compute these vectors at a large 
number of points, and use cubic spline interpolation to 
estimate the rates of change of orbital constants between 
these points. This allows us to compute the temporal 
evolution of all relevant quantities, including the orbital 
trajectories and gravitational waveforms. 

We implemented this algorithm, discrctizing the or- 
bital phase space from an initial separation of r = 
10, 000 mi to the Kerr ISCO 



Hsco 
mi 



3 + Z 2 T Zi)(3 + Z 1 +2Z 2 ), 



Z, = l + (l-g?) 1/3 [(l + ?i) 1/;i + (l-gi) 

Z 2 = {M\ + Z 2 f\ 
in a 1, 000 point grid, equally spaced in 



a/3 



{ mi L0) 1/3 = 



3/2 / 3/2 



4 



J /mf 



1/3 



(22) 



(23) 



We cannot evolve inside the Kerr ISCO with such a 
frequency-domain code, as stable orbits do not exists in 
this regime. We cannot evolve inside the Kerr ISCO with 
a frequency-domain code, as stable orbits do not exist in 
this regime (and so such orbits do not have a well-defined 
frequency spectrum) . The code we use to construct our 
waves is based on [28|, [2^, |6(| , updated to use the spec- 
tral methods introduced by Fujita and Tagoshi (6ll. l62j . 
A detailed presentation of this code and its results is in 
preparation [6^ . 

The dominant error in these Teukolsky-based wave- 
forms is due to truncation of the sums over t and m. 
In particular, to compare with PN results, we must map 
this spheroidal decomposition to a spherical one (see Ap- 
pendix [A] for more details). Such a mapping requires one 
to include a buffer region of £ modes about the largest 
mode computed. We have been careful to use a suffi- 
ciently wide buffer and total number of £ modes such 
that the energy fluxes are accurate to 10~ 10 for all ve- 
locities and spins. In particular, this means that in the 
strong field region (close to the ISCO) up to 50 £ modes 
were included. Other sources of error are due to the 
intrinsic double precision in the numerical solution to 
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the Teukolsky equation, the discretization of the orbital 
phase space, and its cubic-spline interpolation. All of 
these amount to errors of order 10 -14 . All in all and in 
terms of GW phase, the Teukolsky-based waveforms are 
accurate to at least 10~ 2 rads. during an entire year of 
evolution. 



III. CALIBRATING THE TEST-PARTICLE 
ENERGY FLUX 

We consider here a calibrated EOB model that is built 
from the h^ m functions in Eq. (|17|) . but in which higher- 
order PN terms are included in the functions pi m and 
are calibrated to the numerical results. In particular, we 
write 



„22 
PCal 



n 33 



,(9.1) 



i. jo 



622' ^ eulerlog 2 v 2 



(12,0) 



0^2'^ eulerlog 2 ?/ 



„33 



,(8,2) 
^33 



^33 eulerlog 3 ?; 2 



q v 



q 2 v 8 



(10,0) . ,(10,0) 



''33 



'3,3 



eulcrloggU^ 



(24) 



where { a ^ M \b^ M * > ) are eulerlog- independent and 
eulerlog-dependent calibration coefficients that enter the 
(£,m) mode at 0(v N ) and proportional to q M . As in 
Refs. [H, HH| the euler-log function is defined as 



eulerlog m (a;) = j E + log (2my/x) , 



(25) 



where 7^ = 0.577215 ... is Euler's constant. Notice that 
we have introduced 4 calibration parameters in the non- 
spinning sector of the flux and 4 in the spinning sector. 
The spin parameter q denotes here the spin of the back- 
ground. When we neglect mass-ratio terms, we choose 
q = q\. However, when we switch on the mass-ratio 
terms we have an ambiguity on the choice of q. Follow- 
ing Ref. [m we choose q = q, where q is the deformed- 
Kerr spin parameter defined in Appendix [C] Note that 
since q = qi + 0(m2/m{), these choices are identical in 
the test-particle limit, which is when the energy flux is 
calibrated. 

The choice of calibrating function in Eq. (|24|) is rather 
special and requires further discussion. We have cho- 
sen this function so that leading-order corrections in the 
two dominant GW modes, (2, 2) and (3, 3), are included. 
Higher (£, m) modes contribute significantly less to the 
GW and its associated flux. In each mode, we have in- 
cluded the leading-order unknown terms that are both 
q-indcpcndcnt and g-dcpendcnt. Since g-indepcndcnt 
terms in pi m are known to much higher order (5.5PN) 
than the q-dependent ones (4PN) in the test-particle 
limit (24I KM] . l65j |. spin- independent calibration coeffi- 
cients enter at a much higher PN order. The spin- 
dependence of the calibration terms is inferred from 
known terms at lower PN orders. We have investigated 
many functional forms for the calibrating functions, with 



a varying number of degrees of freedom, and found the 
one above to be optimal in the class studied. 

The introduction of 8 additional calibration terms 
might seem like a lot. In Ref. [52j, the knowledge of 
non-spinning terms up to 5.5 PN order was found to be 
crucial to obtain a sufficiently good agreement in the flux. 
Moreover, only 4 additional calibration terms (2 at 6PN 
order in P22 and 2 at 5PN order in P33) were needed to 
reach a phase agreement of 1 rad after two years of evolu- 
tion. Similarly here, we expect that if the remaining 4.5, 
5, and 5.5 PN order terms were calculated in the test- 
particle limit when the central black carries a spin, then 
the flux would also improve, requiring a smaller number 
of calibration parameters. It is quite likely that those 
higher-order PN terms will be computed in the near fu- 
ture, as they involve dramatically less complicated cal- 
culations than PN terms in the comparable-mass case. 

Having in hand an improved GW energy flux carried 
away to infinity, this must be enhanced with expressions 
for the GW energy flux that is absorbed by the back- 
ground BH. We do so here by simply adding the Taylor- 
expanded form of the latter (see Appendix |B|) to the flux 
of Eq. (|15p . The BH absorption terms in the GW energy 
flux depend on polygamma functions, which are compu- 
tationally expensive to evaluate. We have empirically 
found that expanding this function in q <C 1 to 30th 
order is a sufficiently good approximation for our pur- 
poses. When performing computationally expensive cal- 
culations (like the fits described below) we shall employ 
such expansions, but when solving for the orbital phase 
and when computing the waveforms we shall return to 
the full polygamma expressions. 

The resulting EOB energy flux, including BH absorp- 
tion terms, is then calibrated via a two-dimensional, 
least-squares minimization relative to numerical data ob- 
tained from Teukolsky-based calculations. The fitting 
routine is two dimensional because when considering EM- 
RIs about spinning backgrounds, the flux depends on two 
independent variables: the orbital velocity (or frequency 
or separation) and the spin of the background. This, in 
turn, increases the number of points that need to be used 
by more than an order of magnitude to properly calibrate 
Eq. (|24p . In all fits, we have assumed a data variance of 
10 -11 for all velocities and spins and we have required a 
relative accuracy of one part in 10 8 . Since the data is now 
two-dimensional, one must search for a global minimum 
in [q, v) space. After doing so, we find the calibration 
parameters 



(9,1) 
a 22 

,(12,0) 



-3.1092, 
493.08 , 



4 8 3 ' 2) = -17.310, 
Jio,o) 



6 (9,1) 
"22 

(,(12,0) 



'33 



'33 



= -113.01 



,(10,0) 



= -18.786, 
-247.89, 
= 22.500, 
= 28.125, 



(26) 
(27) 
(28) 
(29) 



The computational cost of the calibrations performed 
in this paper is much larger than those carried out in 
Ref. [53| for the following reasons. First, we consider 
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FIG. 1. Fractional difference between PN and Teukolsky- 
based fluxes as a function of velocity for spins q — (0.5, 0.9) 
(top) and q = (—0.9, —0.5,0.0) (bottom). The dotted curves 
employ the Taylor-expanded PN flux with BH absorption 
terms, while the dashed and solid curves use the uncalibrated 
and calibrated p-resummed fluxes with BH absorption terms 
respectively. 



twice as many calibration parameters as in Ref. [52j |. 
increasing the dimensionality of the fitting space. Sec- 
ond, global minimization routines require non-trivial al- 
gorithms that are numerically more expensive than those 
employed in one-dimensional minimizations. Third, the 
amount of data fitted increases by at least one order of 
magnitude, due to the intrinsic bi-dimensionality of the 
problem. Combining all of this, the computational cost 
of performing the calibration is now more than 100 times 
larger than in Ref. [52| . Even then, however, these fits re- 
quire O(10) CPU minutes to complete. Once they have 
been carried out, this calculation does not need to be 
repeated again in the waveform modeling. 

Figure [JJ plots the fractional difference between the 
analytical GW energy flux and that computed with 
Teukolsky-based waveforms as a function of velocity, 
from an initial value of v/c = 0.01 to the veloc- 
ity at the ISCO, for five different spin values: q = 
(—0.9,-0.5,0.0,0.5,0.9). All comparisons are here nor- 
malized to the Newtonian value of the flux -Fkcwt = 
32/5is 2 v 10 . The different curve styles differentiate be- 
tween analytical models: the dotted curves use the total, 
uncalibrated Taylor-expansion; the dashed curves use the 
uncalibrated p-resummcd flux with BH absorption terms; 
the solid curves use the calibrated p-resummcd flux with 
BH absorption terms. Notice that the calibrated model 
does better than the other two by at least two orders of 
magnitude near the ISCO for all spin-values. 

Several interesting conclusions can be drawn from 
Fig. [TJ First, as obtained in Ref. [46| the uncalibrated 
p-resummed model is better than the Taylor-expanded 
version of the flux, by up to nearly an order of magni- 
tude at the ISCO for all spins. In turn, the calibrated 
model is better than the uncalibrated one by one to two 



orders of magnitude near ISCO for all spins. One could 
also calibrate the Taylor-expanded flux (not shown in 
Fig. [TJ , but this would not produce such good agreement 
in the entire (v, q) space. This is clearly because the un- 
calibrated p-resummed model is more accurate than the 
Taylor one, and thus the calibration terms have to do 
less work to improve the agreement. For the calibrated 
Taylor and p-resummed models to become comparable 
in accuracy one would have to include up to at least 16 
calibration coefficients in the Taylor model. 

The inclusion of BH absorption coefficients is crucial 
to obtain good agreement with the full Teukolsky-based 
flux, a result that was not obvious for the case of non- 
spinning EMRIs. Figure [5] plots the fractional differ- 
ence between the uncalibrated EOB GW energy flux and 
Teukolsky-based one as a function of velocity for five dif- 
ferent spin values: q — (—0.9, —0.5, 0.0, 0.5, 0.9), from an 
initial value of v/c = 0.01 at the ISCO, for five differ- 
ent spin values: q = (—0.9,-0.5,0.0,0.5,0.9). For these 
cases, we have wisco = (0.343,0.367,0.408,0.477,0.609). 
The solid curves use the uncalibrated EOB model includ- 
ing the Taylor-expanded BH absorption contributions, 
while the dotted curves do not. For the non-spinning 
case, observe that there is a very small difference (smaller 
than 10 -2 ) between adding the BH absorption terms or 
not. 

For the spinning cases, however, this is not the case. 
For rapidly spinning backgrounds, adding the BH ab- 
sorption terms improves the agreement by an order of 
magnitude. Presumably, resumming these terms in a 
multipolar-factorized manner would improve the agree- 
ment even more. The BH absorption terms play a much 
larger role in the spinning case because spin changes the 
PN order at which absorption enters in the energy flux. 
These terms enter at 4PN order for Schwarzschild black 
holes, but at 2.5PN order for non-zero spin. This change 
of order has a very large and important impact on the 
system's evolution. 

The inclusion of calibration parameters to improve 
the agreement of PN-inspircd fluxes and Teukolsky-based 
ones for EMRIs is certainly not new. In Ref. [HJ , a sim- 
ilar, PN-inspired calibration was carried out for circular- 
inclined orbits (and more generic ones). Before calibra- 
tion, their fluxes were Taylor-expanded to 2PN order and 
included only the contribution that escapes to infinity 
(not the BH absorption terms discussed above). Their 
fit was then done with Teukolsky-data produced by an 
older version of the code used here, which was accurate 
only to one part in 10 6 . Moreover, the fit was done in 
the range r € (5,30)M [v € (0.183,0.436)], so the fitted 
function loses accuracy rapidly outside this regime, par- 
ticularly close to the ISCO. Inside the fitting regime, the 
flux was fitted to an accuracy of 3 x 10 -2 using 45 cali- 
bration coefficients for an inclined, but fixed orbit. The 
accuracy decreases to 0.1 for orbits which get closer to 
the ISCO. 

To fairly compare the results of Ref. [32[ with our re- 
sults which are restricted to circular, equatorial orbits, we 
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FIG. 2. Fractional difference between p-resummed and 
Teukolsky-based fluxes as a function of velocity for spins 
q = (-0.9,-0.5,0.0,0.5,0.9). The dotted curves do not 
include the Taylor-expanded BH absorption contributions, 
while the solid lines do. 

implemented their model and re-calibrated it considering 
only circular, equatorial orbits. Using 45 coefficients, we 
found an accuracy similar to ours at high velocities close 
to the ISCO, but worse at low velocities. This is be- 
cause their fluxes before calibration are not as accurate 
as the one employed here (by including up to 5.5PN order 
terms and BH absorption terms), particularly at low ve- 
locities. It is important to emphasize that by calibrating 
8 parameters instead of 45 we here obtain better flux ac- 
curacies than in Ref. [32j for circular, equatorial EMRIs. 
We could obtain even better accuracy if we were using 
a larger number of calibration coefficients, e.g. using 16 
coefficients the agreement with the Teukolsky-based flux 
would be of O(10- 5 ). 

IV. COMPARISON OF THE GW PHASE AND 
AMPLITUDE 

The comparison of EOB and Teukolsky evolutions re- 
quires that we choose a specific EMRI. We shall here fol- 
low Ref. [HI and choose system parameters that define 
two classes of EMRIs: 

• System I explores a region between orbital separa- 
tions r/M £ (16, 26), which spans orbital velocities 
and GW frequencies in the range v £ (0.2, 0.25) 
and /gw € (0.005,0.01) Hz respectively. Such 
an EMRI has masses mi = 10 5 Mq and m,2 = 
10 Mq for a mass ratio of 10~ 4 and it inspirals for 
~ (6.3-6.7) x 10 5 rads of orbital phase depending 
on the spin. 

• System II explores a region between orbital sep- 
arations r/M £ (H,nsco), which spans orbital 
velocities and GW frequencies in the range v £ 
(0.3,msco) and / GW € (0.001, fg$°) Hz respec- 



tively. Such an EMRI has masses mi = 10 6 Mq 
and 77i2 = 10 Mq for a mass ratio of 10 -5 and it 
inspirals for ~ (1.9-4.5) x 10 5 rads of orbital phase 
depending on the spin. 

The evolution of Sys. I is stopped around an orbital 
separation of 16M, because this coincides with a GW 
frequency of 0.01 Hz, which is close to the end of the 
LISA sensitivity band. The evolution of Sys. II is usu- 
ally stopped at the orbital separation corresponding to 
the ISCO, or whenever its GWs reach a frequency of 
0.01 Hz. For each of these systems, we shall investigate 
different background spin parameters. 

Before proceeding, notice that Sys. I and II should 
not be compared on a one-to-one basis. One might be 
tempted to do so, as Sys. I resembles a weak- field EMRI, 
which inspirals at a larger orbital separation and with 
smaller orbital velocities than Sys. II, a more strong-field 
EMRI. Comparisons are not straightforward, however, as 
these systems accumulate a different total number of GW 
cycles. In fact, Sys. I usually accumulates almost twice 
as many GW cycles as Sys. II. Therefore, even though 
one might expect PN models of Sys. I to agree better 
with Teukolsky-based evolutions, this need not be the 
case, as this system has more time (as measured in GW 
cycles) to accumulate a phase and amplitude difference 
than Sys. II. 

We compare the EOB and the Teukolsky-based wave- 
forms after aligning them in time and phase. Such an 
alignment is done by minimizing the statistic in Eq. (23) 
of Ref. [49j|, just as was done in Ref. (52J. This is equiva- 
lent to maximizing the fitting factor over time and phase 
of coalescence in a matched filtering calculation with 
white noise [49j]. The alignment is done in the low- 
frequency regime, inside the time interval (0, 64)A GW , 
where A GW is the GW wavelength. This quantity de- 
pends on the spin of the background, ranging from 386 M 
(63 M) to 415 M (121 M) for Sys. I (Sys. II). This corre- 
sponds to aligning the initial phase and frequency inside 
a window of length in the range (0.004, 0.01) months de- 
pending on the system and spin of the background. We 
have checked that increasing the size of the alignment 
window does not affect the final phase and amplitude dif- 
ference; for example, for a spin of q = 0.9 and Sys. I, in- 
creasing the alignment window by a factor of two changes 
the final phase difference by 0.002 rads and the relative, 
fractional amplitude agreement by 0.0004%. 

Figure [3] shows the absolute value of the dephasing and 
relative, fractional amplitude difference of the dominant 
(£, m) = (2, 2) mode for both systems and a variety of 
background spins. For Sys. I, the calibrated EOB model 
maintains a 1 radian phase accuracy over at least the 
first 6 months for all spin values, while for Sys. II the 
same phase accuracy is maintained for up to only the 
first 2 months. The amplitude agreement is also excellent 
for all spin values, with better agreement for Sys. I. As 
found in Ref. [52| the GW phase agreement is primarily 
due to the correct modeling of the orbital phase, as the 
former tracks the latter extremely closely; we find that 
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FIG. 3. Absolute value of the dephasing (left) and relative, fractional amplitude difference (right) computed in the calibrated 
EOB model and the Teukolsky-based waveforms for the dominant (£, m) — (2, 2) mode. Different curves correspond to different 
background spin values. 



the difference between the orbital and GW phase over a 
one year evolution is less than 0.1 rads. 

The agreement in the phase as a function of back- 
ground spin follows closely the flux agreement shown in 
Fig. [T] This is hard to see in Fig. [T] which is why Fig. [4] 
zooms into the velocity region sampled by Sys. I and plots 
all background spin cases for the calibrated p-resummed 
system. Observe that the q = 0.0 case has the best flux 
agreement, which explains why the phase and amplitude 
agreement is so good for this case in Fig. [3J Observe 
also that the q = 0.9 and q = —0.9 cases have the worst 
flux agreement, which also explains why they disagree 
the most in phase and amplitude in Fig. [3] 



rads [~ 10 3 -10 4 rads] and - 0.1% (~ 10%) for Sys. I 
(Sys. II) after a one year-evolution for different spin val- 
ues. The above results are consistent with the arguments 
in Ref. [23 |. who concluded that 3.5PN accurate GW 
phase expressions could lead to phase errors around 10 3 - 
10 4 radians over the last year of inspiral. That analysis 
reached those conclusions by comparing 3 to 3.5PN ac- 
curate, analytic expressions for the GW phase. Here, we 
are comparing full-numerical evolutions of the PN equa- 
tions of motion carried out to much higher order, and, 
of course, we find that such conclusions depend sensi- 
tively on the type of EMRI considered and the spin of 
the background. 
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FIG. 4. Fractional difference between the calibrated 
p-resummed and Teukolsky-based fluxes for spins q = 
(-0.9,-0.5,0.0,0.5,0.9) as a function of velocity. We plot 
here only the range of velocities sampled by Sys. I. 

The accuracy of the calibrated EOB model is excellent 
relative to Taylor-expanded PN models. If one were to 
use an uncalibrated Taylor-expanded version of the flux, 
instead of the calibrated p-resummed flux, one would 
find a phase and amplitude disagreement of ~ 10 1 -10 2 



FIG. 5. Absolute value of the dephasing computed for the 
dominant mode in different PN models and the Teukolsky- 
based waveforms. Different curve colors/shades correspond to 
different background spin values, while different curve types 
correspond to different PN models. 

The increase in accuracy of the calibrated p-resummed 
model is due both to the calibration and to the ht m fac- 
torized resummation. This fact can be appreciated in 
Fig. [SJ where we plot the absolute value of the dephasing 
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for the dominant mode in different PN models and the 
Teukolsky-based waveforms. Light curves (orange in the 
color version) correspond to background spins of q = 0.9, 
dark curves (red in the color version) to q = —0.9 and 
black curves to non-spinning backgrounds. Dotted curves 
use the uncalibrated Taylor flux model, dashed curves 
the uncalibrated p-resummed model and solid curves the 
calibrated version. For Sys. I, there is a large gain in ac- 
curacy by switching from the uncalibrated Taylor model 
to the uncalibrated p-resummed model, but then the cal- 
ibration itself does not appear to improve the accuracy 
substantially. For Sys. II, on the other hand, the calibra- 
tion can increase the accuracy up to almost 2 orders of 
magnitude, as in the q = 0.9 case. 

The agreement in the phase and amplitude is not only 
present in the dominant mode, but also in higher (£, m) 
ones, as shown in Fig. [6] We here plot the absolute 
value of the dephasing and the relative, fractional am- 
plitude difference between the calibrated EOB model 
and Teukolsky-based waveforms for the dominant mode, 
as well as the (3,3) and (4,4) modes. We have here 
shifted the higher (£, m) modes using the best frequency 
shift that maximizes the agreement for the dominant 
mode. This agreement is simply a manifestation of the 
agreement in the orbital phase. In fact, we find that 
$gw ~ "^orb, with differences that are always less than 
1 rad for all systems considered. 

Higher-^ modes contribute significantly less to the SNR 
than the dominant (2,2) mode. Figure [7] plots the rela- 
tive fraction between the squared of the SNR computed 
with only the hi m component of the waveform and that 
computed by summing over all modes. This figure uses 
data corresponding to Systems I and II, both with spin 
q = 0.9 (results for other spin values are almost identi- 
cal). Clearly, the (2, 2) mode is dominant, followed by the 
(3,3) and (4,4) modes. Because of this feature of quasi- 
circular inspirals, obtaining agreement for the (2, 2) mode 
implies one can recover over 97% of the SNR. 



V. DATA ANALYSIS IMPLICATIONS 

Although the phase agreement presented in the previ- 
ous section is a good indicator of the validity of the EOB 
model, one is really interested in computing more realistic 
data analysis measures. In this section we compute the 
mismatch between the Tcukolsky and the EOB model, 
maximized over extrinsic parameters and as a function 
of observation time. 

Let us first introduce some basic terminology. Given 
any time series a(t) and b(t), we can define the following 
inner-product in signal space 



(a\ 6) =4 Re 



ajf) b*(f) 

Sn(f) 



(30) 



where the overhead tildes stand for the Fourier transform 
and the star stands for complex conjugation. The quan- 
tity S n (f) is the spectral noise density curve, where here 



we follow [66|, [67J ■ Notice that we use the sky-averaged 
version of this noise curve here, which is larger than the 
non-sky-averaged version by a factor of 20/3. In particu- 
lar, this means that our SNRs are smaller than those one 
would obtain with a non-sky-averaged noise curve by a 
factor of (20/3) 1 / 2 - 2.6. 

Given this inner-product, we can now define some use- 
ful measures. The SNR of signal a is simply 



P = V W a) , 
while the overlap between signals a and b is simply 



(31) 



M 



(«l b) 



\/H «) (b\ b) 



(32) 



with the mismatch MM = 1 — M. The max label here is 
to remind us that this statistic must be maximized over 
a time shift and a phase shift (see eg. Appendix B of [IH 
for a more detailed discussion). 

The data analysis measures introduced above (p and 
MM) depend on the length of the time-series, i.e. the ob- 
servation time. Figure [8] plots the mismatch between the 
Teukolsky-based waveforms and a variety of models for 
both Sys. I and II and a background spin of q = 0.9 as 
a function of observation time. The vertical lines cor- 
respond to observation times of 2 weeks, 2 months, 6 
months, 9 months and 11.5 moths, together with their as- 
sociated SNRs at 1 Gpc. The mismatches are computed 
with different analytical models: black crosses stand for 
the calibrated p-model with 8 calibration coefficients; red 
circles to the uncalibrated p- model; blue squares to the 
uncalibrated Taylor model; green circles to that EOB 
evolution using the original flux of Ref. [HI which has 
45 calibration coefficients (denoted GG in the figure). 
For comparison, we also include the amount of dephas- 
ing (numbers next to data points in Fig. [8]) at 2 weeks, 2 
months, 6 months, 9 months and 11.5 months. Observe 
that the calibrated p-model maintains an overlap higher 
than 97% over 9 and 4 months for Sys. I and II respec- 
tively. The uncalibrated p-model performs comparably 
to the EOB model using the flux of Ref. [32| which has 
45 calibration coefficients, both of which have an over- 
lap higher than 97% over 6 and 1 month for Sys. I and 
II respectively. In the case of Sys. II the calibrated flux 
of Ref. [HJ perform better than the uncalibrated pe m 
model. Also observe that the uncalibrated Taylor model 
is simply inadequate to model EMRIs for any observation 
time. 

We then see that the use of the calibrated EOB model 
allows us to integrate over longer observation times, com- 
pared to a 2-week period or to other models. In turn, this 
allows us to recover a higher SNR that we would other- 
wise. The increase in SNR scales as the square root of 
the observation time, as expected. For example, since 
the calibrated EOB model is accurate over 9 months, 
one would be able to coherently recover an SNR of 64 
at 1 Gpc for Sys. I, to be compared with an SNR of 49 
obtained after 6 months of coherent integration if using 
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t [Months] t [Months] 



FIG. 6. Absolute value of the dephasing (left) and relative, fractional amplitude difference (right) computed in the calibrated 
EOB model and the Teukolsky-based waveforms. The solid curve corresponds to the dominant (2, 2) mode, while the dashed 
curve is for the (3, 3) mode and the dotted curve for the (4, 4) mode. Different curves stand for different background spins. 
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FIG. 7. Relative fraction between the squared of the SNR 
computed with only the (£, m) mode and that computed with 
all {£, m) modes. The lines connecting points are only meant 
to group points with £ — m (dotted), £ = m + 1 (dashed) and 
£ = m + 2 (dot-dashed). 



for example the fluxes of Ref. [32j. In general, an inte- 
gration over a period longer than two weeks gains us a 
large increase in SNR. Such gains in SNR are important 
because they allow us to see EMRIs farther out. Since 
the SNR scales as p ~ D^ 1 , where Dl is the luminosity 
distance, even an SNR increase in a factor of 3 increases 
our accessible volume by a factor of 27, since the later 
scales as D\. 



VI. HIGHER-ORDER EFFECTS 

Let us now discuss how finite mass-ratio higher-order 
effects affect the GW phase and amplitude. Those ef- 
fects are encoded in the v terms present in the radiation- 
reaction force and in the Hamiltonian. The former are 



System 


I 


I 


I 


II 


II 


II 


<7i 


-0.9 


0.0 


0.9 


-0.9 


0.0 


0.9 


No rel. v 


10.04 


1.60 


9.36 


7.63 


42.21 


48.39 


v in H 


30.38 


0.083 


18.96 


4.38 


40.35 


47.13 


v in E 


19.99 


5.31 


4.49 


7.24 


41.50 


46.55 


vinH and E 40.32 


6.83 


14.08 


3.98 


39.64 


45.29 



TABLE I. Absolute value of the total dephasing after a 11.5 
moths of evolution. The first row includes no relative v con- 
tributions in the Hamiltonian or the radiation-reaction force. 
The second row includes relative v terms in the Hamiltonian, 
while the third row includes such terms in the radiation reac- 
tion force. 



second-order effects in the dissipative dynamics, while 
the latter are first-order effects in the conservative dy- 
namics. We have analytic control over the PN version of 
such effects within the EOB formalism, but until now, we 
had set both of these to zero when comparing to Teukol- 
sky evolutions, as the latter do not account for such ef- 
fects. In the EOB model, however, it is straightforward to 
include such terms, as PN expansion are formally known 
for all mass ratios at some given order in v. The inclu- 
sion of conservative v terms is achieved by using the spin 
EOB Hamiltonian of Ref. |4(j reviewed in Appendix \C\ 
The inclusion of dissipative v terms is achieved by in- 
cluding relative v terms in the multipolarly decomposed 
waveform hi m and flux J- . 

Whether such mass-ratio effects matter depends on the 
EMRI considered. In Ref. [52|, it was found that such ef- 
fects increase the dephasing between EOB and Teukolsky 
models by one order of magnitude after a two year evo- 
lution for non-spinning EMRIs, a result consistent with 
Ref. (6^. This effect is greatly amplified when consider- 
ing spinning EMRIs. Table U compares the effect that the 
inclusion of relative v terms in the EOB Hamiltonian and 
the radiation-reaction force has on the final dephasing af- 
ter a one- year evolution. In order to read out the effect 
of such finite mass-ratio terms in the phasing, one must 
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t [Months] t [Months] 

FIG. 8. Mismatch between Teukolsky and EOB waveforms as a function of waveform duration for Systems I and II and a 
background spin of q = 0.9. 



compare rows two, three and four to the baseline given in 
the first row of Table |U For example, the effect of the v 
terms in the Hamiltonian are such as to increase the de- 
phasing by 27.20 — 10.04 = 17.16 radians. Observe that 
the conservative and dissipative v terms usually push the 
dephasing in different directions, partially canceling out 
when both of them are present. Even then though, the 
generic effect of high-order v terms is to increase the rate 
of dephasing by several tens of radians are a one year 
evolution. Notice furthermore that the magnitude of the 
effect is here not very large because we are considering 
circular equatorial orbits. 

Finite mass-ratio effects are clearly suppressed when 
dealing with circular orbits. This is because, for such 
orbits outside the ISCO, the effect of the conservative 
self force is simply to shift the waveform from one or- 
bital frequency to another. Thus, from an observational 
standpoint, such an effect is unmeasurablc. Even though 
the conservative self- force shifts the ISCO, this effect is 
still degenerate with a shift of the system's mass pa- 
rameters. This discussion, however, neglects radiation- 
reaction, which is crucial to model a true inspiral wave- 
form. One can think of the radiation-reaction force as 
defining a trajectory through the sequence of orbital en- 
ergies that an orbit follows. There is gauge invariant 
information in this sequence, in the sense that the map- 
ping between energies and orbital frequencies depends 
on the details of the orbit at each energy level. For a 
given radiation-reaction force, the sequence of geodesic 
orbits (and hence the sequence of frequencies) depends 
on whether the conservative self- force is included or not. 

As is clear from this discussion, the "real" (gauge- 
invariant) effect of the conservative self-force on quasi- 
circular inspiral waveforms can be rather subtle. A ro- 
bust effect, however, does arise if the self-force acts on 
a more generic orbit, such as an eccentric one. In that 
case, this force will act separately on the radial and the 
azimuthal orbital frequencies, which can leave a poten- 
tially strong imprint in the waveform. In principle, even 
for an inclined circular orbit there could be a strong im- 



print. In practice, however, the azimuthal and polar or- 
bital frequencies are quite similar, which suggests that 
perhaps, even in this case, the self-force effects will be 
small. 

With all of this in mind, let us discuss the results pre- 
sented in Table U in more detail. Our study suggests that 
the overall effect of v terms in both H and T leads to only 

5.2 and 2.5 additional radians of phase for non-spinning, 
Sys. I and II respectively. This is in fact consistent with 
the results presented in Ref. (52[, except that there one 
considered 2-ycar long evolutions. One might wonder 
whether using the non-spinning Hamiltonian of Ref. (49| 
(where the deformed-Schwarzschild potential are Padc- 
resummed instead of being given by Eqs. (|C18|) . flC8[) ) 
has an effect on this dephasing. We have investigated 
this question and found that the additional contribution 
to the phase is 0.05 (0.03) and 1.06 (1.33) radians for 
Sys. I and II respectively over the entire year of inspi- 
ral using the 3PN (4PN) Padc form of the deformed 
potentials. [We notice [69| that the 4PN Pade poten- 
tials of Ref. [49[ reproduce very closely the ISCO-shift of 

This implies that the non-spinning Hamilto- 
at 3PN and 4PN order is sufficiently close to the 
Hamiltonian presented in Appendix B for data analysis 
of non-spinning EMRIs. 

One can also compare the results in Table U to the re- 
cent study of Huerta and Gair [7(| , who investigated the 
effect of i/ 2 -corrections in the determination of parame- 
ters, given an EMRI signal. Their Table I presents the 
number of cycles accumulated for a variety of mass ratios. 
Their last column happens to correspond to our Sys. II 
with no spin, for which they get a total dephasing of 

2.3 rads and 3.8 rads after the last year of inspiral when 
including only conservative and all second-order correc- 
tions. This is to be compared to our results: 1.86 rads in 
and 2.6 rads after the last year of inspiral when including 
only conservative and all second-order corrections. These 
numbers are in excellent agreement, allowing for differ- 
ences in the modeling. Their analysis suggests that such 
small difference will not affect parameter estimation for 
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System 


I 


I 


I 


II 


II 


II 




-0.9 


0.0 


0.9 


-0.9 


0.0 


0.9 


q 2 = -0.75 


36.34 


3.96 


16.03 


4.41 


40.41 


47.04 


g 2 = 0.00 


40.32 


6.83 


14.08 


3.98 


39.64 


45.29 


q 2 = 0.75 


44.31 


9.70 


12.14 


3.56 


38.88 


43.54 



TABLE II. Absolute value of the total dephasing in rads after 
a 11.5 moths of evolution. The first row sets the second BH's 
spin to zero, while the second row sets it to q 2 = 0.75. In both 
cases, the spin of the background is set to q\ = 0.9 and v 7^ 
in both the Hamiltonian and the radiation-reaction force. 

EMRIs similar to Sys II. For Sys I, however, the dephas- 
ing is much larger as the mass ratio is less extreme by one 
order of magnitude; thus, parameter estimation might be 
affected in this case. 

Another high-order effect that one can study is the 
inclusion of the small object's spin in the evolution of 
the binary. Since the spin angular momentum of the 
small body scales with its mass, its contribution to the 
orbital evolution is one order in v suppressed. We can 
model this effect by allowing q^ ^ in the EOB Hamil- 
tonian of Ref. [!(| and letting v ^ 0. Doing so for Sys. I 
(Sys. II) with qi = 0.9 and q2 = 0.75, we find that the to- 
tal dephasing now becomes 17.26 (45.49 rads), as shown 
in Table HI1 This is to be compared to the case when 
q-2 = 0, which returns a dephasing of 14.08 rads (45.29 
rads). Thus, the effect of the second spin contributes 
roughly 3.3 rads in this case. We can also compare these 
results to those estimated by Barack and Cutler [l||. In 
their Appendix C, they estimate that for quasi-circular 
inspirals similar to our Sys. II, the spin of the second 
body should induce a dephasing of roughly 1 — 10 rads. 
This is in good agreement with the results corresponding 
to Sys. II in our Table IE Our results are also in good 
agreement with an upcoming and independent investiga- 
tion of spin-effects in the PN phasing |7l[ . 

Finally, taking into account the results of including 
finite mass-ratio effects, we can conclude that unless 
these are precisely modeled, it is not worth requiring 
an agreement better than 10-30 rads when calibrating 
the phase of the test-particle EOB waveforms against the 
Teukolsky-bascd waveforms. 



VII. CONCLUSIONS 

We have constructed an EOB model for EMRIs in 
quasi-circular, equatorial orbits about spinning back- 
grounds. In the test-particle limit, this model consists of 
adiabatically evolving a test-particle in the Ker r sp ace- 
timc using the factorized energy flux of Refs. (45l |46|. 
augmented by 8 calibration coefficients. The latter are 
determined by comparing the factorized energy flux to a 
Teukolsky-based flux, built from solutions to the Teukol- 
sky equation in the radiative approximation. In the adi- 
abatic approximation, the EOB waveforms can be con- 
structed in CPU seconds at a very low computational 



cost. When finite mass-ratio effects and the small ob- 
ject's spin are included, we build the EOB model by 
numerically solving the Hamilton equations with the 
spin EOB Hamiltonian of Ref. [l8j and the Teukolsky- 
calibrated factorized energy flux augmented by finite 
mass-ratio effects [46[ . 

For both EMRI systems considered, we find excellent 
phase and amplitude agreement, with dephasing less than 
one radian, over periods of months. The exact length 
of the agreement depends on how relativistic the EMRI 
system is. We also calculated the overlap between EOB 
and Teukolsky-based waveforms to find it higher than 
97 % over 4 to 9 months, depending on the EMRI system 
considered. 

The EOB waveforms built here have higher overlaps 
and better phase agreements that all currently known 
EMRI models for spinning, equatorial systems, while re- 
quiring much fewer calibration parameters. In particu- 
lar, the EOB model with 8 calibration coefficients out- 
performs by almost an order of magnitude the numerical 
kludge waveforms with the calibrated fluxes of Ref. [32j 
which use 45 calibration coefficients. This implies that 
EOB waveforms with 8 calibration coefficients can be 
used for longer coherent integrations, allowing us to ob- 
tain a 50 % increase in SNR. In turn, this implies that 
EOB waveforms can see EMRIs that arc farther out, in- 
creasing the accessible volume by at least a factor of two 
relative to numerical kludge waveforms. Furthermore if 
we were using 16 calibration coefficients, we could im- 
prove the dephasing from 0.91 rads (8.7 rads) to 0.85 
rads (4.2 rads) for System I (System II) after 6 months 
of evolution. In turn, this would decrease the mismatch 
from 0.2% (12%) to 0.19% (3.9%) for System I and II 
after 6 months of evolution. 

Another possible avenue for future research is the cal- 
culation of high PN order terms in the energy flux and 
waveforms. Our EOB model relies on the use of accurate 
fluxes, but for spinning systems, the flux to infinity is 
only known up to 4PN order in the test-particle limit. 
This is in contrast to the non-spinning terms that are 
known to 5.5PN order or the BH absorption terms that 
are known to 6.5PN order. The calculation of the spin- 
dependent terms in the flux to infinity to 4.5, 5 and 5.5PN 
order terms in the test-particle limit is not quixotic and 
would be invaluable. Once these coefficients are known, 
then presumably the EOB waveforms would be more ac- 
curate and might require less adjustable parameters. 

Of course, the EOB waveforms we constructed here are 
less powerful than kludge waveforms [IH, [H, 113 in their 
generality. Our waveforms cannot yet model inclined or 
eccentric inspirals. The inclusion of inclined orbit should 
be relatively straightforward, but the addition of eccen- 
tricity might require some revamping of the EOB frame- 
work. Future work in this direction would be definitely 
worthwhile. 

Ultimately, one would like to obtain a waveform model 
that is sufficiently fast, efficient and accurate to do real- 
istic EMRI data-analysis for LISA. Such a model would 
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need to include the correct self-force contributions to 
the conservative dynamics, the appropriate second-order 
terms in the radiation-reaction force and the correct 
terms that describe the small object's spin. The EOB 
model we developed here does agree with all known PN 
self-force calculation to date. However, since full self- 
force calculations are not yet completed, we do not have 
a way of assessing the error of including the currently 
known EOB finite mass-ratio effects. In fact, so far the 
only comparisons between the PN/EOB and self- force 
results have been concerned with the non-spinning case, 
and have been limited to the ISCO shift [IH|,[Hl[ and other 
gauge invariant quantities [72|, [73| . Quite interestingly, 
the calibration of the EOB model to comparable-mass 
numerical-relativity simulations improves the agreement 
of the model to self- force results [Hi], |(3{|. Thus, we hope 
that future calibrations of the spin EOB Hamiltonian to 
comparable-mass simulations of spinning BHs will allow 
us to build an EOB model which includes finite mass- 
ratio effects in a more accurate way. 

All that said, we have found that the presence of PN 
self-force terms in the spin EOB model of Ref. [4(| leads 
to dcphasing of 10-30 rad over one year depending on the 
EMRI system and the spin of the background. By con- 
trast, the inclusion of the small object's spin introduces 
dephasing on the order of a few radians. Taking into 
account those results, we can conclude that unless those 
finite mass-ratio effects are precisely modeled, currently, 
it is not worth requiring an agreement better than 10-30 
rads when calibrating the phase of the test-particle EOB 
waveforms against the Teukolsky-based waveforms. This 
in turn implies that calibrating 16 parameters instead of 
8 to an EOB model is overkill as other systematics (in- 
duced by neglecting self-force effects) will be dominant. 
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Appendix A: Mapping between Spheroidal and 
Spherical Harmonics 

1. Quantities at Spatial Infinity 

Let us first consider the GW fluxes that describe radi- 
ation escaping to infinity, and later discuss the radiation 
absorbed by the BH. Frequency-domain Tcukolsky equa- 
tion codes expand the curvature scalar ^4 as 

4>4 = lY, Z "nSi m (8,<t>)e-^ t ■ (Al) 

im 

We have here incorporated the e™^ dependence into the 
spin- weighted spheroidal harmonics ST . These harmon- 
ics depend on the value of migiw m , and the minus super- 
script is a reminder that the we consider here harmonics 
of spin- weight —2. Throughout this appendix, the index 
£ refers to the spheroidal harmonic index, while I refers 
to the spherical harmonic index. 

From ipi, we compute waveforms via ip4 = {h + — 
ih x )/2, and hence, for a frequency-domain application, 

h+-ih* = -^^S-JMK^ 4 • (A2) 

im m 

We have implicitly assumed that u m is time- independent, 
or at least that its time-dependence is subleading. A 
better expansion is to re-express things in terms of the 
accumulated phase, i.e., the integral of the frequency 
w m = $ m , namely 



V'4 = -E Z ^^ n (^,0)e-^ (A3) 

im 




The EOB and numerical relativity (NR) communi- 
ties like to project these quantities onto a basis of 
spin- weighted spherical harmonics. For -04, they define 
Ci m (t,r) via 

^ = ;£3*.(*> r )*h;(M)- (A5) 

Ira 

and the harmonically-decomposed waveforms hi m (t,r) 
via 

h = \Y. h ^ r ) Y i^ 9 ^) ■ ( A6 ) 

Ira 

The minus superscript again denotes spin- weight —2. 
Defining the inner-product 

(y im \f)= J ^Y~ m *(0,^)f , (A7) 

the extraction of the Ci m and hi m is simple: 
C lm (t,r) =r(Y lm \^) , h lm (t,r) =r(Y l ~\h) . (A8) 
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Let us now take the Schwarzschild limit to see whether 
these expressions simplify. In that limit 



57 



= Y: 



(A9) 



and thus, performing the necessary projections and tak- 
ing advantage of the orthonormality of spherical harmon- 
ics, we find 



C lm (t,r) = Z{ 
hi m (t,r) 



9 7 H 



(A10) 
(All) 



In spinning backgrounds, however, the mapping is 
more complicated. We can use the fact that spheroidal 
harmonics can be expressed as a sum of spherical har- 
monics via 



5 



3 



(A12) 



The dependence on miqiuj m enters through the expan- 
sion coefficients bj (see, e.g. Ref. [5!|). Inserting this 
expansion into Eq. (jAlj) and using the inner-product def- 
inition, Eq. (| A8|) . we find 
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b\Z, 
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I 
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(A13) 
(A14) 



In the Schwarzschild limit, = Sfj, so that Kerr simply 
limits as it should. 

From the definition of the Isaacson stress-energy ten- 
sor, one can easily show that 



d 2 E°° _ v 

dtdn ^ ■ 



I 7 H I 2 
^2 If^ml 

47ro; 2 n 



H 12 



Zra V i 



4tuj 2 



(A15) 
(A16) 



We have here used the orthonormality of both the 
spheroidal and spherical harmonics to simplify the sums, 
as well as the fact that the time dependence of the 
disappears when its modulus is computed. Performing 
the angular integrals leaves us with familiar formulas: 
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47TO& 



(A17) 



This breaks down nicely enough that it is useful and sen- 
sible to define the modal contributions Ei m . 



2. Quantities at Event Horizons 

As a general principle, computing quantities that are 
related to an event horizon is usually more complicated 



than computing the same quantities at spatial infinity. 
For the fluxes, for example, this is because there is no 
simple generalization of the Isaacson tensor on the hori- 
zon. Instead, one must examine the shear of the horizon's 
generators, look at how this shear generates entropy, 
and then apply the area theorem to compute fluxes 74 1 . 
The relevant quantities at the horizon depend on the 
Ncwman-Penrose scalar ipo, instead of ip4, the former of 
which is a quantity of spin-weight +2, rather than —2. 

The GW energy flux per unit solid angle at the horizon 
is given by 



d 2 E 



dtdtt 



u m m x r + 1 HH 2 
2ir Pm 1 



(A18) 



where p Tl 



and where uj + = q 1 /{2r + ) is the 



angular velocity of observers co-rotating with the event 
horizon. The quantity cr HH is the shear to the horizon's 
generators as found by Ref. [75|]. This quantity is fairly 
simply related to ipo, so let us introduce an expansion of 
ipo in spin-weight +2 spheroidal harmonics 



-2 \ woo 



Im 



wz(r)st m {e,<j>)> 



(A19) 



where A is given in Eq. ([7]). Notice that this quantity 
diverges on the event horizon r + because the Kinnersley 
tetrad, which is used to define the projection for ipo , is 
ill-behaved as r — > r + . This can be corrected for by 
converting to er HH for any given (£, m) mode [74| 



(A20) 

The complex number "f m is given by j m = — [A{ip m + 
2e)(2mir + ) 2 ] -1 , where e = \J m\ — m 2 g 2 /(4mir + ). 

With this in hand, the GW energy flux at the horizon 
becomes 

J2^^ hmflWiZn st m ) 2 (A21) 



dtdil 



2-npr, 



E ( S im 



\wz\ 2 



tmJ 16p m (p2 i + 4 e2)(2m 1 r + )3 4^ 
which integrates to 



E H = 



E 

im 



I6p m (p m +4e 2 )(2mir + ) 3 4ttw 



W,°°\ 2 

(A22) 



Implementing this equation is difficult because it requires 
that one computes both ipo and ip4 when solving the 
Teukolsky equation. Since these quantities have differ- 
ent angular dependence and a different source function, 
this would be a non-trivial undertaking. 

Instead, one can take advantage of a remarkable sim- 
plification, the so-called Starobinsky identities (76|, that 
relates ip4 to ipo and vice-versa via an algebraic relation. 
We can use this to relate the coefficients to the co- 
efficients Z%L, namely 



(A23) 
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where 



and 



a 6i(2 mi r + )H Pm {p 2 m + 4e 2 )(-i Pm + 4e) voo 

Pirn = 6* 



|Q m | 2 = [(A + 2) 2 + 4q imi uj m - 4q 2 m 2 u 2 m ] (A 2 
+ 3Qm,qiniiUj, n — 36g 2 m 2 w 2 „) + (2A + 3) 
x (96<7 2 m 2 w 2 n — 48mqimiu m ) 
+ lULuf n m 2 (l - q 2 ) , 
Im Q m = l2miio , 



Re C lm = +y\cim\ 2 - IMmjuj 2 , 

A = £i m — 2q\m\muj + q\m\uj 2 — s(s + 1) , (A24) 

and £g m is the eigenvalue of (6,<j)). 

We can then finally write the energy flux formula as 



Ot-lr, 



tin 



(A25) 



where the coefficient ag m , given explicitly in Ref. [591 ] . 
agglomerates the factors /3, 7 into one big expression. 
This expression reorganizes terms slightly in order for its 
structure to resemble the expression for the flux to infin- 
ity as much as possible. Using this, it is simple to write 
the modal contribution as decomposed into a spheroidal 
harmonic basis: Ef m = a £m |Z^ l | 2 /(47ra; 2 a ). 

Decomposing the horizon energy flux formula into a 
spherical harmonic basis is slightly more difficult. The 
key confusing issue is that there are now two spheri- 
cal harmonic basis to worry about: one for each spin- 
weight. With respect to these two bases, we can define 
two spheroidal harmonic expansions: 



(A26) 
(A27) 



(The coefficients dj expand the +2 spin-weight spheroidal 
harmonic in +2 spin-weight spherical harmonics, just as 
the coefficients bj do so for the —2 harmonics.) It's worth 
emphasizing that the different spherical harmonics are 
not simply related to one another. 

The two quantities which can be put into a spheri- 
cal harmonic basis are ■i/'o and tp4, both evaluated in the 
vicinity of the horizon: 



i> a = A- 2 J2WZd e J Y+ m 



= A" 2 V [/°° Y 



(A28) 
(A29) 
(A30) 



1p 4 = 



A 2 



(r — iq\m\ cos#) 4 
A 2 

(r — %q\m\ cos#) 4 



5>^rrJ^)e-fA31) 

ijm 

jm 1 ji 



(A32) 
(A33) 



Using the results presented in this appendix, one can eas- 
ily find an expression for E H in terms of the +2 harmonic 



coefficients U^: 



^I6p m (p 2 m + 4e 2 )(2m ir+ r 4™ 2 



(A34) 



Unfortunately, this is not that useful, as it requires 
knowledge of the ipo expansion coefficients W£fL ■ 

Taking advantage of the Starobinsky identity again, we 
can combine Eqs. (|A23|) and (|A30j) to find 



U°° 



(A35) 



It is then a straightforward to insert this into Eq. ([A"34|) 
to obtain the down-horizon flux expanded into modes of 
+2 spherical harmonics. 



3. Truncation issues 

We have so far been rather schematic regarding the 
limits on all sums. In principle, all these sums should be 
carried out from some lower limit Z m j n to infinity, where 
the former is given by i m ; n = min(|s|, \m\). In a numerical 
application, the upper limit must be truncated at some 
finite value l max . We typically find that the magnitude 
of terms falls off as a power of I. When decomposing 
into spheroidal harmonics, it is thus typically sufficient 
to pick some cutoff value and apply it uniformly. 

Applying such a cutoff is slightly more complicated 
when we convert to spherical harmonics. The reason 
is that a given spheroidal harmonic £ has contributions 
from spherical harmonics at index j > £. Consider, 
as a concrete example, the spheroidal harmonic for 
a/mi = 0.99, u = 0.1: the expansion coefficients for this 
harmonic are 



b% = -0.0110657, b% 



0.99987, 



b\ 



0.0117368, 



0.000123221 , b% = 9.4336 x 10" 



6.4276 x 10" 
1.93317 x 10 



b\ 



3.70511 x 10" 11 , 



-13 



b 5 n 



.97558 x 10 



-16 



(A36) 



Coefficients beyond 6^ are small enough that our code 
does not compute them in this case. Notice that as we 
move away from the j = £ term (whose value is close 
to unity) the coefficients fall off by roughly powers of 
e ~ 0.01. This behavior is typical, although the value 
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of e depends strongly on q\ and u> (e.g., e 
q x = 0.1, w = 0.1, but e ~ 0.4 for q 1 = 0.999, 



10 3 for tions to the radiation flux through the horizon are 



5). 



When one converts from Zf m 



to Cr 



s jm , this behavior 
forces us to include a buffer of £- values beyond the maxi- 
mum spherical harmonic that we want to compute. The 
size of the buffer depends (rather strongly) on the values 
of qi and u. In the weak field, where q\m\uj <C 1 even for 
large to, it is enough to include a buffer of 2 (i.e., such 
that f max = jmax + 2). For the results used in this paper, 
we have chosen a uniform buffer region of size 8. 



Appendix B: GW Energy Absorption 



The iV/2-th order Taylor expansion of the flux in PN 
theory is given by the series [l9( 



N 

4^ = ^Ncwt. M U ) + Iog(«)] V n , (Bl) 

n=0 



where J-Newt = 32/5^ 2 i; 10 is the leading-order (Newto- 
nian) piece of the flux, v is the circular orbital frequency, 
log stands for the natural logarithm while a n and b n are 
PN parameters, with b n <6 = 0. 

The PN parameters can be classified according to their 
physical origin and whether they include spin contribu- 
tions or not. The flux pieces that account for GW emis- 
sion to infinity are well-known and, for example, are given 
in Ref. [zj [Za| ■ Those that correspond to radiation in- 
falling into the horizon will be labeled (a" or ,&^ or ), with 
a superscript S (NS) if they are further spin-dependent 
(spin-independent) . 

The coefficients associated with BH absorption can 
only be formally obtained employing BH perturbation 
theory, as PN theory treats BHs as effective test parti- 
cles. The logarithm-independent terms associated with 
non-spinning contributions to the radiation flux through 
the horizon are 



Hoi.NS -i 

o 8 ' =1, 

Eor.NS i 

a w = 4, 
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(B2) 
(B3) 

(B4) 

(B5) 



where a<°g' N = 0. All logarithm-dependent terms iden- 
tically vanish here: 6" or,NS = 0. Similarly, the logarithm- 
independent terms associated with spinning contribu- 
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- V (0) 
x ^ (0) 
x V (0) 
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-2N-1/2 



3 + (1 - q 

Z-iq{l-q 2 Y 1,2 ])-'Mq{l + ?>q 2 ) 
2\-l/2' 



3 + 2^ (i - r 

3 - 2iq (1 - g 



2\-i/2 



3iq(l + 3f) 

(Bll) 



<5 =0 and where the polygamma function 



where a 

V> (n) ( z ) = (d n T(z)/dz) T^)- 1 is the nth-derivative of 
the Gamma function. The coefficients (a^ r,s , a^'' 5 ) are 
also known, but we do not write them out here as they 
are lengthy and unilluminating [e.g., see Appendix J in 
Ref. 24|]. Notice that the BH absorption coefficients in 
the spinning case are non-zero starting at 2.5 PN order, 
which is to be contrasted with the non-spinning BH ab- 
sorption terms that start at 4 PN order. 

An ambiguity exists when incorporating these BH ab- 
sorption contributions into the flux. As one can ob- 
serve, the spin-dependent coefficients a^ or,s depend on 
the spin parameter of the background q, for which one 
could choose the real spin parameter q = q\ or the effec- 
tive spin parameter q = q, defined in Appendix [U] Since 
q = <7i + 0(m2/mi), these choices are identical in the test 
particle limit, when we calibrate to Teukolsky-fluxes. In 
lack of better guidance, we here choose q = q. 



Appendix C: Spin EOB Hamiltonian 



In Sec. I VII we have investigated how the analyti- 
cal results calibrated to the Teukolsky-based waveforms 
change when we switch on the PN conservative self-force 
and the second-order radiation reaction effects, and when 
wc include the spin of the small object. This study em- 
ployed the spin EOB Hamiltonian of Ref. , which was 
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derived from Ref . (5?| , building also on results of R ef. f39lj . 
As we shall review below, the Hamiltonian of Ref. |4(| re- 
produces the known spin-orbit (spin-spin) PN couplings 
through 2.5PN (2PN) order for comparable masses, and 
all PN couplings linear in the spin of the small object in 
the test-particle limit. We shall here restrict attention to 
circular, equatorial orbits, and assume that the spins are 
aligned with the orbital angular momentum. 

The motion of a spinning test-particle in a generic 
curved spacetime is described by the Papapetrou equa- 
tion froMsiT ] . Reference [13] derived a Hamiltonian whose 
Hamilton equations are equivalent to the Papapetrou 
equation. This Hamiltonian therefore describes the mo- 
tion of a spinning particle in a generic curved spacetime, 
and Ref. [4(| computed it in the particular case of the 
Kerr spacetime in Boycr-Lindquist coordinates. Denot- 
ing with Si and mi the spin and mass of the background 
BH, and with S2 and m.2 the spin and the mass of the 
smaller BH, the Hamiltonian of a spinning test-particle 
in Kerr has the generic form floj 



H = Hi 



NS 



H s , 



(CI) 



where J/ns is the Hamiltonian of a non-spinning test par- 
ticle in Kerr, given by Eq. ([5]), while Hg depends on Si 



and S2 and, if PN expanded, generates all PN terms lin- 
ear in the small object's spin S2. 

In Ref. |4(j the authors constructed the spin EOB 
Hamiltonian by mapping the PN Hamiltonian of two 
BHs of masses 7711,2 and spins Si^ into the effective 
Hamiltonian of a spinning test-particle of mass /1 = 
m,\mij {rri\ + 1112) and spin S* moving in a deformed- 
Kerr spacetime with mass M = mi + TTI2 and spin Skcit, 
v = n/M being the deformation parameter. Note that 
the deformed-Kerr spin parameter q = |SKerr|/Af 2 7^ q\, 
but instead q ps qi(l — 2m2/mi + . . .) when m2jvn\ <C 1. 

The effective Hamiltonian is 1401 



NS 



/' 



2M r : 



(C2) 



where -Hns is the Hamiltonian of a non-spinning effective 
particle in the deformed-Kerr background, 



(C3) 



which differs from Eq. (|9]) in that the Kerr potentials 
A and Wfd have been replaced by their deformed forms 
At and Wfd (and also mi — > M, qi — > q, rri2 — > /it). 
Furthermore, Hg in Eq. (|C2[) is linearly proportional to 
the effective particle's spin S* and reads 




where we denote with a comma the derivative with re- 
spect to r. The term proportional to S 2 in Eq. (|C2[) is 
added to reproduce known spin-spin results at 2PN or- 
der. The quantities (|C3[) and (|C4[) depend on the Kerr- 
deformcd potentials A t , A r , At, uj[d, while 



Q = l 



H 2 A t 



In particular, we have [401 ] 

At = (r 2 + M 2 q 2 ) 2 - M 2 q 2 A t , 

and 

n »r2 fd fd 9 3 ^ 4 

w fd = 2q M r + u\ d v - h v 

r r 



(C5) 



(C6) 



(C7) 



where Lo id and Lo f 2 A are two adjustable parameters regulat- 
ing the frame dragging strength. Although precise values 



for u>i and uj 2 can only be determined by calibrating 
the model against NR simulations of comparable-mass 
spinning BHs, a preliminary comparison of the final spin 
predicted by the EOB model to NR results [H H sug- 
gests that io[ d ps -10 and iof ps 20. 

The deformed-Kerr potential A t is given at 3PN order 

by 



At = r 2 [A(u) + q 2 u 2 ] 



A(u) = 1 -2u + 2vu 3 + v ( — - — 7T 2 ) u 4 

3 32 



(C8) 



(C9) 



When setting v = 0, A t reduces to the Kerr expression 
(with mi -> M, gi g) A = A f = r 2 -2M r+g 2 7\/ 2 . 
In order to guarantee the presence of deformed horizons 
(which correspond to the zeros of A t ), Ref. suggested 
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to rewrite Eq. (fUS|) as (u = M/r) 
'l-2u(l- Ku) 



, 2 2 

+ g u 



with 



and 



A 
Ai 



A 



A, 



[\-Kvf 
1 + i/Ao+z/log (l + J2 A 

X = A' + 4(Aq - A ) i/ , 



(CIO) 
(Cll) 

(C12) 



A' A - 2) , 
-2(^A-1)(A + A ), 

i Ai (-4i/ A + Ax + 4) - q 2 (v A - l) 2 A 



(C13) 
(C14) 



3 — 77 



- A\ + 3(v K — 1) A 2 + 3A 2 Ax 
6(vK -l)(-vK + A 2 + l) 
3q 2 ( l /A-l) 2 A 1 l , 



(C15) 



(C16) 



^{6, 2 (A 2 -2A 2 )(,A-1) ; 

-8(i/A-l)A?-12A 2 A 2 
+12[2(^A-1)A 2 + A 3 ] Ax 

+6 [A 2 -4A 3 (i/A'-l)] }, 



3A? 



(C17) 



When expanding Eq. (|C10|) through 3PN order, one re- 
covers Eq. (|C8[) . The quantity (|C10|) depends on two 
parameters Ao and Ax- Ao is fixed to the value 1.4467 
in order to reproduce the results of Ref. [25[ for the shift 
of the ISCO frequency due to the conservative part of the 
self force. Also, recent comparisons of the EOB model 
with numerical simulations of non-spinning comparable 
mass BHs have suggested Ai sa 3/4. 



The deformed-Kerr potential A r is given by (4CJ 

A r = A t { 1 + log[l + 6u u 2 + 2(26 - 3v) vu 3 ]} , 

(C18) 

which reduces to the Kerr-potential A in the limit v = 
(with mi — > M. qi — > g). 

Finally, the spins Sjcerr and S* in the effective descrip- 
tion are not equal to Si a 5^2, are instead given by 



S* = S 1 ^ + S 2 ^ + ±A S , 
mi Tri2 c 

<JKerr = S\ + S2 , 



(C19) 
(C20) 



where 



12 l r \ 771 1 7Jl2 



-4(Si + S 2 ) 



3(Si + 5 2 ) 



-4(Sx^ + S 2 ^ 
mi m 2 



(C21) 



With all of this at hand, the EOB Hamiltonian used in 
Sec. EE is 



H 



EOB 



= M 




(C22) 



A few final observations are due at this point. When 
the smaller BH has zero spin S* 2 = and mass m 2 <C 
mi, at lowest order in m 2 /mi the EOB Hamiltonian of 
Eq. (|C22j) reduces to the Hamiltonian of a non-spinning 
test particle in Kerr. This is because both S* and the de- 
formations of the Kerr potentials are C(m 2 /mi). How- 
ever, at the next-to-lcading order in the mass-ratio, the 
EOB Hamiltonian presents corrections with respect to 
the Hamiltonian of a non-spinning test particle in Kerr, 
(i) because of the deformations of the Kerr potentials; (ii) 
because of the effective spin S*, which is not zero; (Hi) 
because of the higher-order terms in v that one obtains 
expanding Eq. (|C22|) . These corrections encode the con- 
servative part of the self- force in the EOB framework ■ 
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